#Nonlinear SDOF time history analysis
#DEVELOPED BY CM

wipe

set duration CSZ;

#location loop
for {set location 1} {$location < 10.5} {incr location} {
	
	set timestepfile [open "C:/Program Files (x86)/OpenSees/Column/$duration/GM/$location/timestep.txt" "r"];                  
	set timestep [read $timestepfile];
	close $timestepfile;

	set gmlengthfile [open "C:/Program Files (x86)/OpenSees/Column/$duration/GM/$location/gm_length.txt" "r"];                  
	set gmlength [read $gmlengthfile];
	close $gmlengthfile;
		
	#ground motion loop
	for {set gm 1} {$gm < 60.5} {incr gm} {
		puts "gm $gm"
			
		set dt [lindex $timestep [expr {$gm-1}]];
		set steps1 [lindex $gmlength [expr {$gm-1}]];
		set steps2 [expr {$steps1+30/$dt}];
		set steps [expr int($steps2)];
			
		#period loop
		for {set i 1} {$i < 25.5} {incr i} {
			if {$i==1} {set T 0.10};
			if {$i==2} {set T 0.20};
			if {$i==3} {set T 0.30};
			if {$i==4} {set T 0.40};
			if {$i==5} {set T 0.50};
			if {$i==6} {set T 0.60};
			if {$i==7} {set T 0.70};
			if {$i==8} {set T 0.80};
			if {$i==9} {set T 0.90};
			if {$i==10} {set T 1.00};
			if {$i==11} {set T 1.10};
			if {$i==12} {set T 1.20};
			if {$i==13} {set T 1.30};
			if {$i==14} {set T 1.40};
			if {$i==15} {set T 1.50};
			if {$i==16} {set T 1.75};
			if {$i==17} {set T 2.00};
			if {$i==18} {set T 2.25};
			if {$i==19} {set T 2.50};
			if {$i==20} {set T 2.75};
			if {$i==21} {set T 3.00};
			if {$i==22} {set T 3.25};
			if {$i==23} {set T 3.50};
			if {$i==24} {set T 3.75};
			if {$i==25} {set T 4.00};
				
			set g 386.4;											# gravity
			set m 1;												# mass
			set k		[expr {$m/$T/$T*4*3.14*3.14}];				# stiffness
			set kp		[expr 0.05*$k];								# post-yield stiffness
			set Fy		[expr {max(0.07,min(0.2/$T,0.86))*$m*$g}];	# yield strength
			set negFy	[expr {-1*$Fy}];							# yield strength (negative)
			set Dy		[expr {$Fy/$k}];							# yield displacement
			set Dp		[expr {5*$Dy}];								# plastic displacement at peak strength
			set Fpeak	[expr {$Fy+$kp*$Dp}];						# peak strength
			set Fr		[expr {0.01*$Fy}];							# residual strength
			set Dpc		[expr {($Fpeak-$Fr)/(0.1*$k)}];				# post-peak deformation capacity
			set Dr		[expr {$Dy+$Dp+$Dpc}];						# deformation at residual strength
			set Du		[expr {100*$Dr}];							# ultimate deformation
			set lambda	[expr {100*$Dy}];							# deterioation factor

			# model dimensions and degrees of freedom
			model BasicBuilder -ndm 1 -ndf 1

			# nodes
			node 1 0 
			node 2 0 

			# mass
			mass 2 $m

			# boundary conditions
			fix 1 1
					
			# IMK model:
			uniaxialMaterial Bilin 1 $k 0.05 0.05 $Fy $negFy $lambda $lambda $lambda $lambda 1.0 1.0 1.0 1.0 $Dp $Dp $Dpc $Dpc 0.01 0.01 $Du $Du 1.0 1.0

			#elements
			element zeroLength 1 1 2 -mat 1 -dir 1

			#recorders:
			recorder Node -file IMK/$duration/SDOF/location($location)_T($T)_gm($gm)_disp.txt -node 2 -dof 1 disp
			recorder Node -file IMK/$duration/SDOF/location($location)_T($T)_gm($gm)_force.txt -node 1 -dof 1 reaction
			
			record

			# analysis -------------------------------------------------------------
			# create load pattern
			timeSeries Path 1 -dt $dt -filePath Column/$duration/GM/$location/gm($gm).th -factor [expr ($g)]; # define acceleration vector from file (dt=0.005 is associated with the input file gm)
			pattern UniformExcitation 1 1 -accel 1;		         # define where and how (pattern tag, dof) acceleration is applied
												
			# set damping based on first eigen mode
			set freq [expr [eigen -fullGenLapack 1]**0.5]
			set period [expr {2*3.14/$freq}]
			#puts "frequency is $freq"
			#puts "period is $period"
			set dampRatio 0.05
			rayleigh [expr 2*$dampRatio*$freq] 0. 0. 0.

			# create the analysis
			wipeAnalysis
			constraints Plain;     				 # how it handles boundary conditions
			numberer Plain;					     # renumber dof's to minimize band-width (optimization), if you want to
			test NormDispIncr 1e-6 20;			 # test
			system BandGeneral;					 # how to store and solve the system of equations in the analysis
			algorithm Newton;					 # use Newton algorithm for linear analysis
			integrator Newmark 0.5 0.25;	     # determine the next time step for an analysis
			analysis Transient;					 # define type of analysis: time-dependent
			analyze $steps $dt;					 # apply 11000 0.01-sec time steps in analysis; includes 30 seconds of free vibration after end of gm
				
			#wipe the analysis
			wipeAnalysis
			wipe

			puts "analysis complete"
		}
	puts "gm $gm"
	}
puts "location $location"
}
